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Abstract 

It is crucially important to investigate effects of temperature on magnetic properties such as 
critical phenomena, nucleation, pinning, domain wall motion, coercivity, etc. The Landau-Lifshitz- 
Gilbert (LLG) equation has been applied extensively to study dynamics of magnetic properties. 
Approaches of Langevin noises have been developed to introduce the temperature effect into the 
LLG equation. To have the thermal equilibrium state (canonical distribution) as the steady state, 
the system parameters must satisfy some condition known as the fluctuation-dissipation relation. 
In inhomogeneous magnetic systems in which spin magnitudes are different at sites, the condition 
requires that the ratio between the amplitude of the random noise and the damping parameter 
depends on the magnitude of the magnetic moment at each site. Focused on inhomogeneous mag¬ 
netic systems, we systematically showed agreement between the stationary state of the stochastic 
LLG equation and the corresponding equilibrium state obtained by Monte Carlo simulations in 
various magnetic systems including dipole-dipole interactions. We demonstrated how violations of 
the condition result in deviations from the true equilibrium state. We also studied the characteris¬ 
tic features of the dynamics depending on the choice of the parameter set. All the parameter sets 
satisfying the condition realize the same stationary state (equilibrium state). In contrast, different 
choices of parameter set cause seriously different relaxation processes. We show two relaxation 
types, i.e., magnetization reversals with uniform rotation and with nucleation. 

PACS numbers: 75.78.-n 05.10.Gg 75.10.Hk 75.60.Ej 
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I. INTRODUCTION 


The Landau-Lifshitz-Gilbert (LLG) equation 1 has been widely used in the study of dy¬ 
namical properties of magnetic systems, especially in micromagnetics. It contains a relax¬ 
ation mechanism by a phenomenological longitudinal damping term. The Landau-Lifshitz- 
Bloch (LLB) equation^ contains, besides the longitudinal damping, a phenomenological 
transverse damping and the temperature dependence of the magnetic moment are taken 
into account with the aid of the mean-field approximation. Those equations work well in 
the region of saturated magnetization at low temperatures. 

Thermal effects are very important to study properties of magnets, e.g., the amount of 
spontaneous magnetization, hysteresis nature, relaxation dynamics, and the coercive force in 
permanent magnets. Therefore, how to control temperature in the LLG and LLB equations 
has been studied extensively. To introduce temperature in equations of motion, a coupling 
with a thermal reservoir is required. For dynamics of particle systems which is naturally 
expressed by the canonical conjugated variables, i.e., ( q,p ), molecular dynamics is performed 
with a Nose-Hoover (NH) type reservoir 3 5 or a Langevin type reservoir 6 . However, in the 
case of systems of magnetic moments, in which dynamics of angular momenta is studied, NH 
type reservoirs are hardly used due to complexity 7 . On the other hand, the Langevin type 
reservoirs have been rather naturally applied 2 8 18 although multiplicative nois^® requires the 
numerical integration of equations depending on the interpretation, i.e., Ito or Stratonovich 
type. 

To introduce temperature into a LLG approach by a Langevin noise, a fluctuation- 
dissipation relation is used, where the temperature is proportional to the ratio between 
the strength of the fluctuation (amplitude of noise) and the damping parameter of the 
LLG equation. For magnetic systems consisting of uniform magnetic moments, the ratio is 
uniquely given at a temperature and it has been often employed to study dynamical prop¬ 
erties, e.g., trajectories of magnetic moments of nano-particles 8 } relaxation dynamics in a 
spin-glass system 20 or in a semiconductor 21 . The realization of the equilibrium state by 
stochastic LLG approaches by numerical simulations is an important issue, and it has been 
confirmed in some cases of the Heisenberg model for uniform magnetic moments. ^ 2 ® 
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In general cases, however, magnetic moments in atomic scale have various magnitudes of 
spins. This inhomogeneity of magnetization is important to understand the mechanisms of 
nucleation or pinning . 24 28 To control the temperature of such systems, the ratio between the 
amplitude of noise and the damping parameter depends on the magnetic moment at each 
site. In order to make clear the condition for the realization of the canonical distribution 
as the stationary state in inhomogeneous magnetic systems, we review the guideline of the 
derivation of the condition in the Fokker-Planck equation formalism in the Appendix A. 

Such a generalization of the LLG equation with a stochastic noise was performed to study 
properties of the alloy magnet GdFeCcP^, in which two kinds of moments exist. They ex¬ 
ploited a formula for the noise amplitude, which is equivalent to the formula of our condition 
A (see Sec 0 - They found surprisingly good agreements of the results between the stochas¬ 
tic LLG equation and a mean-field approximation. However, the properties in the true 
canonical distribution is generally different from those obtained by the mean-field analysis. 

The LLG and LLB equations have been often applied for continuous magnetic systems or 
assemblies of block spins in the aim of simulation of bulk systems, but such treatment of the 
bulk magnets tend to overestimate the Curie temperature 1 -^, and it is still under develop¬ 
ment to obtain properly magnetization curves in the whole temperature region 2 11 — 18 . The 
influence of coarse graining of block spin systems on the thermal properties is a significant 
theme, which should be clarified in the future. To avoid such a difficulty, we adopt a lattice 
model, in which the magnitude of the moment is given at each magnetic site. 

Within the condition there is some freedom of the choice of parameter set. In the present 
paper, in particular, we investigate the following two cases of parameter sets, i.e., case A, 
in which the LLG damping constant is the same in all the sites and the amplitude of the 
noise depends on the magnitude of the magnetic moment at each site, and case B, in which 
the amplitude of the noise is the same in all the sites and the damping constant depends on 
the magnitude of the moment, (see Sec 0 )- We confirm the realization of the equilibrium 
state, i.e., the canonical distribution in various magnetic systems including critical region by 
comparison of magnetizations obtained by the LLG stochastic approach with those obtained 
by standard Monte Carlo simulations, not by the mean-held analysis. We study systems 
with not only short range interactions but also dipole-dipole interactions, which causes 
the demagnetizing held statically. We fold that different choices of the parameter set which 
satishes the fluctuation-dissipation relation give the same stationary state (equilibrium state) 
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even near the critical temperature. We also demonstrate that deviations from the relation 
cause systematic and significant deviations of the results. 

In contrast to the static properties, we find that different choices of parameter set cause 
serious difference in the dynamics of the relaxation. In particular, in the rotation type 
relaxation in isotropic spin systems, we find that the dependences of the relaxation time on 
the temperature in cases A and B show opposite correlations as well as the dependences of the 
relaxation time on the magnitude of the magnetic moment. That is, the relaxation time of 
magnetization reversal under an unfavorable external field is shorter at a higher temperature 
in case A, while it is longer in case B. On the other hand, the relaxation time is longer for 
a larger magnetic moment in case A, while it is shorter in case B. We also investigate the 
relaxation of anisotropic spin systems and find that the metastability strongly affects the 
relaxation at low temperatures in both cases. The system relaxes to the equilibrium state 
from the metastable state by the nucleation type of dynamics. The relaxation time to the 
metastable state and the decay time of the metastable state are affected by the choice of 
the parameter set. 

The outline of this paper is as follows. The model and the method in this study are ex¬ 
plained in Sec [TTJ Magnetization processes as a function of temperature in uniform magnetic 


systems are studied in Sec III Magnetizations as a function of temperature for inhomoge¬ 


neous magnetic systems are investigated in Sec. IV, in which not only exchange interactions 
(short-range) but also dipole interactions (long-range) are taken into account. In Sec. |V| 
dynamical aspects with the choice of the parameter set are considered, and the dependences 
of the relaxation process on the temperature and on the magnitude of magnetic moments 


are also discussed. The relaxation dynamics via a metastable state is studied in Sec. VI 


Sec. VII is devoted to summary and discussion. In Appendix [A] the Fokker-Planck equation 
for inhomogeneous magnetic systems is given both in Stratonovich and Ito interpretations, 
and Appendix [B] presents the numerical integration scheme in this study. 


II. MODEL AND METHOD 


As a microscopic spin model, the following Hamiltonian is adopted, 




(ti) 


i^k 


ik 


r 2 
1 ik 
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Here we only consider a spin angular momentum S t for a magnetic moment Mi at each 
site (i is the site index) and regard Mi = S t ignoring the difference of the sign between 
them and setting a unit: g/i-Q = 1 for simplicity, where g is the g-factor and /x B is the Bohr 
magneton 30 . Interaction J t j between the ith and jth magnetic sites indicates an exchange 
coupling, (i,j) denotes a nearest neighbor pair, is an anisotropy constant for the ith 
site, hi is a magnetic field applied to the ith site, and the final term gives dipole interactions 
between the ith and kth sites whose distance is ry*,, where C = is defined using the 
permeability of vacuum /i 0 . 

The magnitude of the moment M % is dehned as Mj = M t |, which is not necessarily 
uniform but may vary from site to site. In general, the damping parameter may also have 
site dependence, i.e., and thus the LLG equation at the ith site is given by 


d 

dt 


Mi 


*7 Mi x Hf + 


cq „ , 

—Mi x 

Mi 


dMi 

dt 


( 2 ) 


or in an equivalent formula: 


d 

dt 


Mi 


7 

1 + a? 


Mi x Hf 


(1 + afMi 


Mi x (Mi x Hf), 


( 3 ) 


where 7 is the gyromagnetic constant. Here Hf is the effective field at the ith site and 
described by 

d 

Hf = , M N ,t ] (4) 

, which contains fields from the exchange and the dipole interactions, the anisotropy, and 
the external field. 

We introduce a Langevin-noise formalism for the thermal effect. There have been several 
ways for the formulation to introduce a stochastic term into the LLG equation. The stochas¬ 
tic field can be introduced into the precession term and/or damping ternl 8 9 11 . Furthermore, 
an additional noise term may be introduced 1012 . In the present study we add the random 
noise to the effective field Hf —» Hf + ^ and we have 


d 

dt 


Mi 


7 

1 + a 2 


MiX(Hf +eo 


Oii 7 


(1 + a 2 ) Mi 


Mi x (Mi x (Hf + &)), 


( 5 ) 


where f is the fi(= 1,2 or 3 for x,y or z) component of the white Gaussian noise applied at 
the ith site and the following properties are assumed: 


(m) = 0 , (e k m(f) = 2 - *)• 


( 6 ) 
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We call Eq. ^ stochastic LLG equation. We derive a Fokker-Planck equation 6 8 for 
the stochastic equation of motion in Eq. (|5j) in Stratonovich interpretation, as given in 
appendix [A] 


i 


7 d f 

1 + a? ' \ 


-7 DiMi x (Mi x 


a 


JjMi x (Mi x Hf) 


d 


dM; 


P(M\, • • • , M n , t) > . 


( 7 ) 


Here we demand that the distribution function at the stationary state (t —>■ oo) of the 
equation of motion (Eq. (JT]) ) agrees with the canonical distribution of the system (Eq. (JTj) ) 
at temperature T, i.e., 


P eq (M 1 , • • • , M n ) oc exp ^ /3'H(M l , • • • , M N ) ), (8) 

where (3 = 

Considering the relation 

- ,M N )=f3HfP eq (M 1 ,--- , M n ), (9) 

we find that if the following relation 


a 


rrjr ~ 7A/3 = 0 
Mi 


( 10 ) 


is satisfied at each site i, the canonical distribution in the equilibrium state is assured. 

When the magnetic moments are uniform, i.e., the magnitude of each magnetic moment 
is the same and M* = \M, t = M, the parameters cq and Di are also uniform cq = a and 


Di = D for a given T. However, when M* are different at sites, the relation (10) must be 
satisfied at each site independently. There are several ways of the choice of the parameters 
oti and Di to satisfy this relation. Here we consider the following two cases: A and B. 

A: we take the damping parameter cq to be the same at all sites, i.e., ai — a 2 — ■ ■ ■ — = 

a. In this case the amplitude of the random held at the ith site should be 


a k B T 1 
Dj =-oc — 

Mi 7 Mi 


( 11 ) 


B: we take the amplitude of the random held to be the same at all sites, i.e., A = A = 
■ ■ ■ — Dm = D. In this case the damping parameter at the ith site should be 

D'yMi 


Oti = 


knT 


oc Mi. 


( 12 ) 


7 










1 


0.8 

0.6 
E 

0.4 

0.2 

0 

0 1 2 3 4 5 6 

T 

FIG. 1: (color online) Comparison of the temperature dependence of m in the stationary state 
between the stochastic LLG method and the Langevin function (green circles). Crosses and boxes 
denote m in case A (ct = 0.05) and case B (D = 1.0), respectively. In the stochastic LLG 
simulation A t = 0.005 was set and 80000 time steps (40,000 steps for equilibration and 40,000 
steps for measurement) were employed. The system size N = L 3 = 10 3 was adopted. 

We study whether the canonical distribution is realized in both cases by comparing data 
obtained by the stochastic LLG method with the exact results or with corresponding data 
obtained by Monte Carlo simulations. We set the parameters 7 = 1 and k-g = 1 hereafter. 

III. REALIZATION OF THE THERMAL EQUILIBRIUM STATE IN HOMOGE¬ 
NEOUS MAGNETIC SYSTEMS 

A. Non-interacting magnetic moments 

As a first step, we check the temperature effect in the simplest case of non-interacting 
uniform magnetic moments, i.e., Jij = 0, Vf- — 0, C — 0 in Eq. ([lj and M* = M (or 
Si = S ), where a and D have no site f-dependence. In this case the magnetization in a 
magnetic held (h) at a temperature (T) is given by the Langevin function: 



* 

m 

m 

- H 

m 

m 

m 

m 



We compare the stationary state obtained by the stochastic LLG method and Eq. (13). 






We investigate m(T ) at h = 2 for M = 1. Figure [I] shows m{T ) when a = 0.05 is fixed (case 
A) and when £) = 1.0 is fixed (case B). We find a good agreement between the results of 
the stochastic LLG method and the Langevin function in the whole temperature region as 
long as the relation (10) is satisfied. Numerical integration scheme is given in Appendix [B| 
The time step of A t = 0.005 and total 80000 time steps (40000 steps for equilibration and 
40000 steps for measurement) were adopted. 


B. Homogeneous magnetic moments with exchange interactions 


Next, we investigate homogenous magnetic moments (M l 
mensions. The following Hamiltonian ((7 = 0, Jij = J, T) A 

fl): 


\M,\ = M ) in three di- 
V A , and h(t) = h in Eq. 


(14) 


h = ~y, js ‘- s i - E ® A ( S ? ) 2 - E hS > 

{id) i i 

is adopted. 

There is no exact formula for magnetization (m) as a function of temperature for this 
system, and thus a Monte Carlo (MC) method is applied to obtain reference magnetization 
curves for the canonical distribution because MC methods have been established to obtain 
finite temperature properties for this kind of systems in the equilibrium state. Here we 
employ a MC method with the Metropolis algorithm to obtain the temperature dependence 
of magnetization. 

In order to check the validity of our MC procedure, we investigated magnetization 
curves as functions of temperature (not shown) with system-size dependence for the three¬ 


dimensional classical Heisenberg model (V A = 0 and h = 0 in Eq. (14)), and confirmed that 


the critical temperature agreed with past studies 3 ^, where /crT c = 1.443J for the infinite 
system size with M — 1. 

We give m{T ) for a system of M — 2 with the parameters J — 1, h — 2 and T> A = 1.0 
for cases A and B in Fig [ 2 } The system size was set N = L 3 = 10 3 and periodic boundary 
conditions (PBC) were used. Green circles denote m obtained by the Monte Carlo method. 
At each temperature (T) 10,000 MC steps (MCS) were applied for the equilibration and 
following 10,000—50,000 MCS were used for measurement to obtain m. Crosses and boxes 
denote m in the stationary state of the stochastic LLG equation in case A (a = 0.05) and 
in case B (D = 1.0), respectively. Here At = 0.005 was set and 80000 steps (40000 for 
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FIG. 2: (color online) Comparison of temperature (T) dependence of m between the Monte Carlo 
method (green circles) and the stochastic LLG method in the homogeneous magnetic system with 
M = 2. Crosses and boxes denote case A with a = 0.05 and case B with D = 1.0, respectively. 

transient and 40000 for measurement) were used to obtain the stationary state of m. The 
m{T) curves show good agreement between the MC method and the stochastic LLG method 
in both cases. We checked that the choice of the initial state for the MC and the stochastic 
LLG method does not affect the results. The dynamics of the stochastic LLG method leads 
to the equilibrium state at temperature T. 


IV. REALIZATION OF THE THERMAL EQUILIBRIUM STATE IN INHOMO¬ 
GENEOUS MAGNETIC SYSTEMS 


A. Inhomogeneous magnetic moments with exchange interactions 


Here we study a system which consists of two kinds of magnitudes of magnetic moments. 


The Hamiltonian (14) is adopted but the moment M, = \Mi\ has i-dependence. We investi¬ 


gate a simple cubic lattice composed of alternating M — 2 and M — 1 planes (see Fig. [3] (a)), 
where J = 1, h = 2 and D A = 1.0 are applied. We consider two cases A and B mentioned 
in Sec. HH 

The reference of m(T ) curve was obtained by the MC method and is given by green 
circles in Figs. [3] (b) and (c). In the simulation, at each temperature (T) 10,000 MGS were 
applied for the equilibration and following 10,000—50,000 MGS were used for measurement. 
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FIG. 3: (color online) (a) A part of the system composed of alternating M = 2 (red long 

arrows) and M = 1 (short blue arrows) layers, (b) Comparison of temperature (T) dependence of 
m between the Monte Carlo method (green circles) and the stochastic LLG method for a = 0.05. 
At = 0.005 and 80,000 steps (40,000 for transient time and 40,000 for measurement) were employed. 
Crosses denote m when Di = D(Mi ) = was used. Triangles and Diamonds are m for 

Di = D( 1) = — for all i and D, = D( 2) = for all i, respectively, (c) Comparison 

of temperature (T) dependence of m between the Monte Carlo method (green circles) and the 
stochastic LLG method for D = 1.0. At = 0.005 and 80,000 steps (40,000 for transient time and 
40,000 for measurement) were employed. Crosses denote m when a* = a(Mj) = was used. 

Triangles are m for ai = a(Mj = 1) = for all i and Diamonds are m for at = a(2) = D ^ b j 2 

for all i. 
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The system size N = L 3 = 10 3 was adopted with PBC. In case A, a(= 0.05) is common for 
all magnetic moments in the stochastic LLG method and Mj (or S)) dependence is imposed 
on Dj as D, = D(Mj) = In case B, D — 1.0 is common for all magnetic moments in 

the stochastic LLG method and cq = a(Mj) = ■ Crosses in Figs. |3] (b) and (c) denote 

m by the stochastic LLG method for cases A and B, respectively. For those simulations 
At = 0.005 and 80,000 steps (40,000 for transient time and 40,000 for measurement) were 
employed at each temperature. In both Figs. [3] (b) and (c), we find good agreement between 
m{T) by the stochastic LLG method (crosses) and m{T) by the MC method (green circles). 

Next, we investigate how the results change if we take wrong choices of parameters. We 
study m(T) when a uniform value Di = D for case A (cq = a for case B) is used for all 
spins, i.e., for both Mi = 1 and M i: = 2. If D(Mi = 2) = |is used for all spins, m{T ) 
is shown by Diamonds in Fig. |3j (b), while if D(M t = 1) = a is applied for all spins, 
m(T ) is given by triangles in Fig. [3] (b). In the same way, we study m(T ) for a uniform 
value of a. In Fig. [3] (c) triangles and diamonds denote m(T ) when a* = a(Mj = 1) and 
a, = a(Mi = 2) are used, respectively. We find serious difference in m(T ) when we do not 
use correct Mj-dependent choices of the parameters. The locations of triangle (diamond) at 
each temperature T are the same in Figs. [3] (a) and (b), which indicates that if the ratio 
a/D is the same in different choices, the same steady state is realized although this state is 
not the true equilibrium state for the inhomogeneous magnetic system. Thus we conclude 
that to use proper relations of M r dependence of Di or cq is important for m{T ) curves of 
inhomogeneous magnetic systems and wrong choices cause significant deviations. 


B. Critical behavior of Inhomogeneous magnetic moments 


In this subsection, we examine properties near the critical temperature. Here we adopt 
the case of h = 0 and V A = 0 in the same type of lattice with M = 1 and 2 as Sec. IV A We 


investigate both cases of the temperature control (A and B). The Hamiltonian here has 0(3) 
symmetry and m is not a suitable order parameter. Thus we define the following quantity 
as the order parameter 31 : 


= \/ml + ml + ml, 


(15) 
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FIG. 4: (color online) Comparison of temperature (T) dependence of m a between the MC method 
(green circles) and the stochastic LLG method for the system of inhomogeneous magnetic moments. 
N = L 3 = 20 3 . PBC were used. In the MC method 10,000 MGS and following 50,000 MGS were 
used for equilibration and measurement at each temperature, respectively. The stochastic LLG 
method was performed in case A with a = 0.05 (erases) and in case B with D = 1.0 (diamonds). 
Here At = 0.005 was applied and 240,000 steps were used (40,000 for transient and 200,000 for 
measurement). 


where 


TV TV TV 

m x = m V=Jf(E, S i), and m z= m = 

2=1 2=1 2=1 


(16) 


In Fig. |4j green circles denote temperature (T) dependence of m a given by the MC 
method. The system size N = L 3 = 20 3 with PBC was adopted and in MC simulations 
10,000 MCS and following 50,000 MCS were employed for equilibration and measurement, 
respectively at each temperature. The magnetizations of m a obtained by the stochastic LLG 
method for case A (crosses) and case B (diamonds) are given in Fig. [4j Here a = 0.05 and 
D = 1.0 were used for (a) and (b), respectively. At = 0.005 was set and 240,000 steps 
(40,000 for transient and 200,000 for measurement) were applied. 

In both cases m a (T) curve given by the stochastic LLG method shows good agreement 


with that obtained by the MC method. Thus, we conclude that as long as the relation (10) 
is satisfied, the temperature dependence of the magnetization is reproduced very accurately 
even around the Curie temperature, regardless of the choice of the parameter set. 
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FIG. 5: (color online) Comparison of temperature (T) dependence of m between the Monte Carlo 
method (green circles) and the stochastic LLG method. Crosses and diamonds denote case A 
with a = 0.05 and case B with D = 1.0, respectively. A reduction of m from fully saturated 
magnetization is observed at around T = 0 due to the dipole interactions. As a reference, m by 
the MC method without the dipole interactions (C = 0) is given by open circles. 

C. Inhomogeneous magnetic moments with exchange and dipole interactions 


We also study thermal effects in a system with dipole interactions. We use the same 
lattice as in the previous subsections. The system is (Jjj = J, V A = V A , and h t (t) = h in 
Eq. Q) given by 


n 


E JS • ■ s- 

(*J> 


E V\SlY - E hS‘ + E^( S -' S ‘- 

■ / 1 ' ih X 

i i i^k ^ 


3 (jPih * * S k ) 


. (17) 


Here a cubic lattice with open boundary conditions (OBC) is used. Since J is much larger 
than C/a 3 (J C/a 3 ) for ferromagnets, where a is a lattice constant between magnetic 
sites. However, we enlarge dipole interaction as C = 0.2 with a = 1 for J = 1 to highlight 
the effect of the noise on dipole interactions. We set other parameters as h = 0.1, V A = 0.1. 
Studies with realistic situations will be given separately. 

We study cases A (a = 0.05) and B (D — 1.0) for this system. We depict in Fig. [5] 
the temperature (T) dependences of m with comparison between the MC (green circles) and 
stochastic LLG methods. Crosses and diamonds denote m(T) for cases A and B, respectively. 
Dipole interactions are long-range interactions and we need longer equilibration steps, and 
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we investigate only a small system with N = L 3 = 6 3 . In the MC method 200,000 MCS 
were used for equilibration and 600,000 steps were used for measurement of m, and for 
the stochastic LLG method At = 0.005 was set and 960,000 steps (160,000 and 800,000 
time steps for equilibration and measurement, respectively) were consumed. A reduction of 
m from fully saturated magnetization is observed. As a reference, m by the MC method 
without the dipole interactions (C = 0) is given by open circles in Fig. [5j This reduction of 
m is caused by the dipole interactions. 

We find that even when dipole interactions are taken into account in inhomogeneous 
magnetic moments, suitable choices of the parameter set leads to the equilibrium state. 
Finally, we comment on the comparison between the LLG method and the Monte Carlo 
method. To obtain equilibrium properties of spin systems, the Monte Carlo method is more 
efficient and powerful in terms of computational cost. It is much faster than the stochastic 
LLG method to obtain the equilibrium m(T ) curves, etc. For example, it needs more than 
10 times of CPU time of the MC method to obtain the data for Fig. [5} However, the MC 
method has little information on the dynamics and the stochastic LLG method is used to 
obtain dynamical properties because it is based on an equation of motion of spins. Thus, it 
is important to clarify the nature of stochastic LLG methods including the static properties. 
For static properties, as we saw above, the choice of the parameter set, e.g., cases A and 
B, did not give difference. However, the choice gives significant difference in dynamical 
properties, which is studied in the following sections. 


V. DEPENDENCE OF DYNAMICS ON THE CHOICE OF THE PARAMETER 
SET IN ISOTROPIC SPIN SYSTEMS (P A =0) 


Now we study the dependence of dynamics on the choice of parameter set. The temper¬ 
ature is given by 


k n T = 


lD t M t 


OLi 


(18) 


which should be the same for all the sites. In general, if the parameter D (amplitude of 
the noise) is large, the system is strongly disturbed, while if the parameter a (damping 
parameter) is large, the system tends to relax fast. Therefore, even if the temperature is 
the same, the dynamics changes with the values of D and a. When the anisotropy term 


exists, i.e., T> A ^ 0, in homogeneous systems (M t = M ) given by Eq. (14), the Stoner- 
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Wohlfarth critical field is h c = 2 MT> A at T = 0. If the temperature is low enough, the 
metastable nature appears in relaxation. On the other hand, if T is rather high or V A = 0, 
the metastable nature is not observed. In this section we focus on dynamics of isotropic spin 
systems, i.e., V A =0. 


A. Relaxation with temperature dependence 


In this subsection we investigate the temperature dependence of magnetization relaxation 
in cases A and B. We adopt a homogeneous system (M* — M — 2) with V A = 0 in Eq. (14). 
Initially all spins are in the spin down state and they relax under a unfavorable external 
held h = 2. The parameter set M = 2, a = 0.05, D = 0.05 gives T = 2 by the condition 


(Eq. (10)). Here we study the system at T = 0.2,1,2, and 10. We set a = 0.05 in case A 
and the control of the temperature is performed by D, i.e. D = 0.005, 0.025, 0.05, and 0.25, 
respectively. In case B we set D = 0.05, and the control of the temperature is realized by 
a, i.e., a = 0.5,0.1,0.05, and 0.01, respectively. 

We depict the temperature dependence of m{t) for cases A and B in Figs. [6] (a) and (b), 
respectively. Here the same random number sequence was used for each relaxation curve. 
Red dash dotted line, blue dotted line, green solid line, and black dashed line denote T = 0.2, 
T = 1, T = 2 and T = 10, respectively. Relaxation curves in initial short time are given in 


the insets. 

In case A, as the temperature is raised, the initial relaxation speed of m becomes faster 
and the relaxation time to the equilibrium state also becomes shorter. This dependence is 
ascribed to the strength of the noise with the dependence D oc T, and a noise with a larger 
amplitude disturbs more the precession of each moment, which causes faster relaxation. 

On the other hand, in case B, the relaxation time to the equilibrium state is longer at 
higher temperatures although the temperature dependence of the initial relaxation speed of 
m is similar to the case A. In the initial relaxation process all the magnetic moments are 
in spin-down state (Sf ~ —2). There the direction of the local held at each site is given 
by ~ Sj + ^ = — 2 x 6 + 2 = —10, which is downward and the damping term 
tends to fix moments to this direction. Thus, a large value of the damping parameter at a 
low temperature T (a oc suppresses the change of the direction of each moment and the 
initial relaxation speed is smaller. However, in the relaxation process thermal fluctuation 
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FIG. 6: (color online) (a) Time dependence of the magnetization (m(t)) in case A, where a = 0.05 
for a homogeneous system with M = 2. Red dash dotted line, blue dotted line, green solid line, 
and black dashed line denote T = 0.2, T = 1, T = 2 and T = 10, respectively. Inset shows the time 
dependence of m(t) in the initial relaxation process, (b) Time dependence of the magnetization 
( m(t )) in case B, where D = 0.05 for a homogeneous system with M = 2. Correspondence between 
lines and temperatures is the same as (a). 


causes a deviation of the local field and then a rotation of magnetic moments from — z 


to z direction advances (see also Fig. 11 ). Once the rotation begins, the large damping 
parameter accelerates the relaxation and finally the relaxation time is shorter. 


B. Relaxation with spin-magnitude dependence 


Next we study the dependence of relaxation on the magnitude of magnetic moments 
in cases A and B. Here we adopt a homogeneous system (Mj = M) without anisotropy( 
T> A = 0) at T = 2 and h = 2. The initial spin configuration is the same as the previous 
subsection. Because 


D oc 


T 

M’ 


and 


M 


a oc —. 

T ’ 


(19) 


raising the value of M is equivalent to lowering temperature in both cases A and B and it 
causes suppression of relaxation in case A, while it leads to acceleration of relaxation in case 
B. Because M affects the local held from the exchange energy at each site, changing the 
value of M under a constant external held h is not the same as changing T and it may show 
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FIG. 7: (color online) Comparison of the time dependence of m between cases A and B by the 
stochastic LLG method. Red and blue lines denote cases A and B, respectively, (a) a = 0.05 for 
case A and D = 1.0 for case B, (b) a = 0.2 for case A and D = 1.0 for case B. 


some modified features. 

In the relation (19), T = 0.2, 1, 2, 10 at M = 2 (Fig|6] (a) and (b)) are the same as 
M = 20, 4, 2, 0.4 at T — 2, respectively. We studied the relaxation ratio defined as m(t)/M 
with M dependence at T = 2 for these four values of M, and compared with the relaxation 
curves of Figj6] (a) and (b). We found qualitatively the same tendency between relaxation 
curves with M dependence and those with 1/T dependence in both cases. A difference was 
found in the initial relaxation speed (not shown). When M > 2, the initial relaxation at 
T — 2 is slower than that of the corresponding T at M = 2. The downward initial local 
field at each site is stronger for larger M due to a stronger exchange coupling, which also 
assist the suppression of the initial relaxation. 

It is found that the relaxation time under a constant external hied becomes longer as 
the value of M is raised in case A, while it becomes shorter in case B. This suggests that 
different choices of the parameter set lead to serious difference in the relaxation dynamics 
with M dependence. 
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VI. DEPENDENCE OF DYNAMICS ON THE CHOICE OF THE PARAMETER 
SET IN ANISOTROPIC SPIN SYSTEMS ( V A / 0) 


A. Different relaxation paths to the equilibrium in magnetic inhomgeneity 


If the anisotropy term exists T> K ^ 0 but the temperature is relatively high, metastable 
nature is not observed in relaxation. We consider the relaxation dynamics when M, has 
i dependence in this case. We study the system (alternating M = 2 and M = 1 planes) 
treated in Sec. |IV A[ We set a configuration of all spins down as the initial state and observe 


relaxation of m in cases A and B. In Sec. IV A we studied cases A (a=0.05) and B (£>=1.0) 


for the equilibrium state and the equilibrium magnetization is m ~ 0.95 at T = 5. We 
give comparison of the time dependence of m between the two cases in Fig. [7] (a), with the 
use of the same random number sequence. The red and blue curves denote cases A and 
B, respectively. We find a big difference in the relaxation time of m and features of the 
relaxation between the two cases. 

The parameter values of a and D are not so close between the two cases at this tempera¬ 
ture (T = 5), i.e., D(M = 1) = 0.25 and D(M = 2) = 0.125 for case A and a(M = 1) = 0.2 
and ct(M = 2) = 0.4 for case B. Thus, to study if there is a difference of dynamics even 
in close parameter values of a and D between cases A and B at T = 5, we adopt common 
a = 0.2, where D(M = 1) = 1 and D(M = 2) = 0.5, as case A and common D = 1.0, where 
a(M = 1) = 0.2 and a(M = 2) = 0.4, as case B. We checked that this case A also gives the 
equilibrium state. In Fig. [T] (b), the time dependence of m for both cases is given. The red 
and blue curves denote cases A and B, respectively. There is also a difference (almost twice) 
of the relaxation time of m between cases A and B. Thus, even in close parameter region of 
a and D , dynamical properties vary depending on the choice of the parameters. 


B. Relaxation with nucleation mechanism 

In this subsection we study a system with metastability. We adopt a homogeneous 
system (M = 2) with J = 1, V A = 1 and h = 2. Here the Stoner-Wohlfarth critical field 
is h c = 2M'Z> V =4, and if the temperature is low enough, the system has a metastable state 
under h = 2. 

At a high temperature, e.g., T = 10 (a = 0.05, D = 0.25), the magnetization relaxes 
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FIG. 8: (color online) (a) Dashed line shows m(t ) for a = 0.05, D = 0.25, and T = 10. Blue 
and green solid lines give m(t) for a = 0.05 at T = 3.5 (case A) and D = 0.25 at T = 3.5 (case 
B), respectively. These two lines were obtained by taking average over 20 trials with different 
random number sequences. The 20 relaxation curves for cases A and B are given in (b) and (c), 
respectively. 


without being trapped as depicted in Fig[8|a) with a black dotted line. When the tempera¬ 
ture is lowered, the magnetization is trapped at a metastable state. We observe relaxations 
in cases A and B, where a = 0.05 for case A and D = 0.25 for case B are used. In Figs, [sjdj) 
and (c), we show 20 samples (with different random number sequences) of relaxation pro¬ 
cesses at T — 3.5 for case A (a = 0.05, D = 0.0875) and case B (D — 0.25, a = 0.143), 
respectively. The average lines of the 20 samples are depicted in Fig |8](a) by blue and green 
solid lines for cases A and B, respectively. In both cases, magnetizations are trapped at a 
metastable state with the same value of m (m ~ —1.55). This means that the metastabil¬ 
ity is independent of the choice of parameter set. Relaxation from the metastable state to 
the equilibrium is the so-called stochastic process and the relaxation time distributes. The 
relaxation time in case A is longer. If the temperature is further lowered, the escape time 
from the metastable state becomes longer. In Figs. [9] (a) and (b), we show 20 samples of 
relaxation at T = 3.1 for cases A and B, respectively. There we find the metastable state 
more clearly. 

Here we investigate the initial relaxation to the metastable state at a relatively low 


temperature. In Figs. 10 (a) and (b), we depict the initial short time relaxation of 20 
samples at T = 2 in cases A (a = 0.05, D = 0.05) and B (D = 0.25, a = 0.25), respectively. 
The insets show the time dependence of the magnetization in the whole measurement time. 
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FIG. 9: (a) and (b) illustrate 20 relaxation curves for a = 0.05 at T = 3.1 (case A) and D = 0.25 
at T = 3.1 (case B), respectively. Metastability becomes stronger than T = 3.5. No relaxation 
occurs in all 20 trials in (a), while five relaxations take place in 20 trials in (b). 


We find that the relaxation is again faster in case B. 

The metastability also depends on M as well as T> A and large M gives a strong metastabil¬ 
ity. Here we conclude that regardless of the choice of the parameter set, as the temperature 
is lowered, the relaxation time becomes longer due to the stronger metastability, in which 
larger D (larger a) gives faster relaxation from the initial to the metastable state and faster 
decay from the metastable state. 

Finally we show typical configurations in the relaxation process. When the anisotropy 
V A is zero or weak, the magnetization relaxation occurs with uniform rotation from — z 
to z direction, while when the anisotropy is strong, the magnetization reversal starts by a 


nucleation and inhomogeneous configurations appear with domain wall motion. In Figs. 11 


we give an example of the magnetization reversal of (a) the uniform rotation type (magneti¬ 
zation reversal for V A = 0 with D = 0.05, T = 2, a = 0.1, M — 4) and of (b) the nucleation 
type (magnetization reversal for V A = 1 with D = 0.25, T = 3.1, a = 0.161, M = 2). 


VII. SUMMARY AND DISCUSSION 

We studied the realization of the canonical distribution in magnetic systems with the 
short-range (exchange) and long-range (dipole) interactions, anisotropy terms, and magnetic 
fields by the Langevin method of the LLG equation. Especially we investigated in detail the 
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FIG. 10: Initial relaxation curves of magnetization. Insets show m(t) in the whole measurement 
time, (a) and (b) illustrate 20 relaxation curves for a = 0.05 at T = 2 (case A) and D = 0.25 at 
T = 2 (case B), respectively. 
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FIG. 11: (a) Typical uniform rotation type relaxation observed in the isotropic spin system, (b) 
Typical nucleation type relaxation observed in the anisotropic spin system. 


thermal equilibration of inhomogeneous magnetic systems. We pointed out that the spin- 
magnitude dependent ratio between the strength of the random field and the coefficient of the 
damping term must be adequately chosen for all magnetic moments satisfying the condition 


(10). We compared the stationary state obtained by the present Langevin method of the 
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LLG equation with the equilibrium state obtained by the standard Monte Carlo simulation 
for given temperatures. There are several choices for the parameter set, e.g., A and B. We 
found that as long as the parameters are suitably chosen, the equilibrium state is realized as 
the stationary state of the stochastic LLG method regardless of the choice of the parameter 
set, and the temperature dependence of the magnetization is accurately produced in the 
whole region, including the region around the Curie temperature. 

We also studied dynamical properties which depend on the choice of the parameters. We 
showed that the choice of the parameter values seriously affects the relaxation process to 
the equilibrium state. In the rotation type relaxation in isotropic spin systems under an 
unfavorable external field, the dependences of the relaxation time on the temperature in 
cases A and B exhibited opposite correlations as well as the dependences of the relaxation 
time on the magnitude of the magnetic moment. The strength of the local field in the initial 
state strongly affects the speed of the initial relaxation in both cases. 

We also found that even if close parameter values are chosen in different parameter sets 
for inhomogeneous magnetic systems, these parameter sets cause a significant difference of 
relaxation time to the equilibrium state. In the nucleation type relaxation, the metastability, 
which depends on V A and M, strongly affects the relaxation in both cases A and B. Lowering 
temperature reinforces the metastability of the system and causes slower relaxation. The 
relaxation to the metastable state and the decay to the metastable state are affected by the 
choice of the parameter set, in which larger D causes fast relaxation at a fixed T. 

In this study we adopted two cases, i.e., A and B in the choice of the parameter set. 
Generally more complicated dependence of M, or T on the parameters is considered. How 
to chose the parameter set is related to the quest for the origin of these parameters. It 
is very important for clarification of relaxation dynamics but also for realization of a high 
speed and a low power consumption, which is required to development of magnetic devices. 
Studies of the origin of a have been intensively performed 32 11 . To control magnetization 
relaxation at finite temperatures, investigations of the origin of D as well as a will become 
more and more important. We hope that the present work gives some useful insight into 
studies of spin dynamics and encourages discussions for future developments in this field. 
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Appendix A: Fokker-Planck equation 


The LLG equation with a Langevin noise (Eq. (j5j)) is rewritten in the following form for 
/i component (fJ. = 1, 2 or 3 for x, y or z ) of the ith magnetic moment, 


dMf 


dt 




Here /■ and g f are given by 


ft = 


7 


1 + a • 


+ ^ X 6 Xpa M"M?Hf’° 


and 


Ll\ 

9i = 


7 


- rT ^ ^AMr + H (-M^ + Mf M A ) 


(Al) 


(A2) 


(A3) 


where H° ,A can have an explicit time (£) dependence, and e Mi ,A denotes the Levi-Civita 
symbol. We employ the Einstein summation convention for Greek indices (/x, zv - - - ). 

We consider the distribution function F = F(AT 1 , • • • ,M^,t) in the 3iV-dimensional 
phase space (M/, Mf , , • • • , M^, M^, Mff ). The distribution function F(iVf 1 , • • • , M N , t ) 

satisfies the continuity equation of the distribution: 

N rl ( y 1 

(A4) 


Ft 




Substituting the relation (Al), the following differential equation for the distribution func¬ 
tion F is obtained. 

N 

rtf. „ a 

(AS) 


a 




2=1 


Regarding the stochastic equation (Al) as the Stratonovich interpretation, making use 


of the stochastic Liouville approach^, and taking average for the noise statistics (Eq. 
we have a Fokker-Planck equation. 

N 


d_ 

Ft 


r) ( 

i =1 * '• 


a/3 ^ / cr/3 


DM? 


(9i p ) 


(A6) 


where P = P(-ZVf 1 , • • • , -Mjv, t) is the averaged distribution function (F). 
Substituting the relation 


d 


dM .! 


<7/3 

:0i = 


7a * -4Aff 


Mi(l + a?) 


(AT) 
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and Eq. (A3) into gf (), we find 


aft / 


d 


'dM\ 


a 9?)= 0. 


(AS) 


Thus Eq.(A6) is simplified to 


d_ 

di 


r\ f r\ 

P(Mi, ,M N ,t) = ~Y,^{(f?~D i gfgf m ')Py (A9) 


Substituting Eqs. (A2) and (A3), we have a formula in the vector representation. 
d 


— P(M \, • • • , M Nl t ) — 


(A10) 


\ ^ 7 9 

it j 


1 + a? dM, 


x (Mj x 


M i X iff + x (Mi x iff) 


5 


3Mi' 


- ,Mtv, 2) f . 


Since ■ (Mj x iff) = 0, it is written as 


9 P(M 1 ,---,Mj V ,t)=^ 7 


dt 


1 + at 

2 L 

-7 DiMi x (Mj x 


a. 


JM, x (Mj x iff) 


(AH) 


d 


dM,: 


P(M 1; • • • , M/v, t) 


In the case that Eq. (Al) is given under Ito definition, we need Ito-Stratonovich trans¬ 
formation, and the corresponding equation of motion in Stratonovich interpretation is 


dM.f 


d 9 i V {Mi) 


~= ,M N ,t)-D l g^(M i r^;;r J +gr(M i )tm- (A 12 ) 


di dM x 

Then the Fokker-Planck equation in Ito interpretation is 

dg\ 


N 




dt 


2—1 


dM? 


d 


(f a - D a Xu — _ P>q a/ q°^ - )P ( 

v* dM x 11 dM a ’ ' 


,\vdgr __ 

dM? 


Since ^ BA/fX 1+^? 


M", the vector representation is given by 


- p ( M 1 ''"’ M - t )=Err^a4 


dt 


O'i 


jj-M, x (M; x H 


eff \ 


- 2 qPjMj - 7AM x (Mj x 


d 


dM j 


P(M 1 ,-- - ,Mw,f)j . 

(A13) 
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Appendix B: Numerical integration for stochastic differential equations 


In stochastic differential equations, we have to be careful to treat the indifferentiability 
of the white noise. In the present paper we regard the stochastic equation, e.g., Eq. (|5]), as 
a stochastic differential equation in Stratonovich interpretation: 

dM? = - + + + (Bl) 

where dW//(t) = J t+dt dst/ 1 / (s) , which is the Wiener process. This equation is expressed by 

dMj 1 — ft*(Mi, ■ ■ ■ ,M N ,t)dt + g? v (M i (t))oc£Wy(t), (B2) 

where o indicates the usage of the Stratonovich definition. 

A simple predictor-corrector method called the Heun method®^, superior to the Euler 
method, is given by 

M?(t + At) = M?(t) 

+ ,M N (t+ At), t + At) + ff(Mi ((),■•• 

+ At)) + g r(M,(t))]AW?, (B3) 

where AW" = W”(t + At) — W(t) and M?(t + At) is chosen in the Euler scheme: 

K(t + At) = M?(t) + ftmit), • • • , M N (t),t)At + g r(Mi(t))AWr (B4) 

This scheme assures an approximation accuracy up to the second order of AW and At. Sev¬ 
eral numerical difference method^ 9 for higher-order approximation, which are often compli¬ 
cated, have been proposed. 

Here we adopt a kind of middle point method equivalent to the Heun method. 

Mf (t + At) = M?(t) 

+ fi(Mt\(t + At/2 ), • • • , M N (t + At/2 ), t + At/2)At 
+ g?'(M i (t + At/2))AWy, (B5) 

where M/‘{t. + At/2) is chosen in the Euler scheme: 

M?(t + At/2) = M?(t) + f?(M 1 (t),--- ,M N (t),t)At/2 + g?(Mi(t))AW", (B6) 
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where A w” = W?{t + At/2) — W"(t). Considering the following relations, 


{AW* AW?) = (\W?{t + At/ 2 ) - wy{t)]\w?(t + At) - wy{t)]) = A At, (B7) 

(A W”) = 0 and (A W”) = 0, this method is found equivalent to the Heun method. We can 


formally replace Aid/ by AWT/2 in Eq. (B6) in numerical simulations. 
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